Mode coupling of Schwarzschild perturbations: Ringdown frequencies 
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Within linearized perturbation theory, black holes decay to their final stationary state through 
the well-known spectrum of quasinormal modes. Here we numerically study whether nonlineari- 
ties change this picture. For that purpose we study the ringdown frequencies of gauge-invariant 
second-order gravitational perturbations induced by self-coupling of linearized perturbations of 
Schwarzschild black holes. We do so through high-accuracy simulations in the time domain of 
first and second-order Regge-Wheeler-Zerilli type equations, for a variety of initial data sets. We 
consider first-order even-parity {I = 2, m — ±2) perturbations and odd-parity (£ — 2, m = 0) ones, 
and all the multipoles that they generate through self-coupling. For all of them and all the initial 
data sets considered we find that — in contrast to previous predictions in the literature — the numer- 
ical decay frequencies of second-order perturbations are the same ones of linearized theory, and we 
explain the observed behavior. This would indicate, in particular, that when modeling or searching 
for ringdown gravitational waves, appropriately including the standard quasinormal modes already 
takes into account nonlinear effects. 

PACS numbers: 04.25.dk, 04.40.Dg, 04.30.Db, 95.30.Sf 
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I. OUTLOOK AND MOTIVATION 

Black hole no-hair theorems [l|, [|| state that within 
Einstein's theory the end point of any system with 
enough gravitational energy to form a black hole is re- 
markably simple: it is uniquely characterized by one 
member of the Kerr family 1 Q, which is described by 
only two parameters: the spin and mass of the final black 
hole. 

As a consequence, the details by which different sys- 
tems decay to such endpoints have been of interest for 
many decades. Pioneering studies were done by a number 
of authors in the early seventies, starting with studies of 
linearized perturbations of non-rotating (Schwarzschild) 
black holes (e.g., 041] )■ Press realized that there is al- 
ways an intermediate stage where the ringdown is dom- 
inated by a set of oscillating and exponentially decaying 
solutions, quasinormal modes (QNMs), whose spectrum 
depends only on the mass of the black hole and the multi- 
pole index I of the initial perturbation . This regime is 
followed by a power-law 'tail' decay due to backscattering 

a- 

In the case of gravitational perturbations of non- 
rotating black holes the relevant equations from which 
QNMs can be inferred are the Regge- Wheeler Q and 



Zerilli [10|, ones. For rotating black holes the corre- 
sponding one (though based on a curvature formalism, as 
opposed to a metric one) is the Tcukolsky equation . 
Their QNMs were first studied by Teukolsky and Press 
[lH . See [HI, [H| for comprehensive reviews on the rich 
area of QNMs. 

The QNM with the lowest frequency is called the fun- 
damental one. Since the subsequent ones (overtones) de- 
cay much faster, the ringdown of Kerr black holes in lin- 
earized theory is in practice described by a few oscillating 
modes which decay exponentially in time, till they reach 
the tail regime. It is interesting to note that the tail 
decay problem for rotating black holes is still not com- 
pletely understood 

From an observational point of view this universal ring- 
down spectrum is of great power: one can use a sin- 
gle QNM detection to infer the mass and spin of the 
black hole source, assuming General Relativity to be cor- 
rect. Alternatively, through a two- mode detection one 
can test General Relativity and/or the assumption that a 
black hole is the source of the measured signal 2l|. The 



1 Charge is expected not to play a significant role in most astro- 
physical scenarios. 



main idea is that the QNM frequencies of both detec- 
tions have to be consistent with respect to their inferred 
masses/spins. 

The LISA mission is expected to measure gravitational 
waves in the low-frequency spectrum: (10~ 5 — 10 _1 ) 
Hz, such as those emitted in the collision of supermas- 
sive binary black holes (SMBBHs) [22|]. Flanagan and 
Hughes [23| showed that, quite generically, the signal to 
noise ratio for these sources in the inspiral regime should 
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be comparable to that one in the ringdown. Therefore, 
detection of SMBBHs by LISA through the measurement 
of QNMs seems to be feasible. Assuming a lower cutoff 
of (1(T 4 - 1CT 5 ) Hz and requiring that the QNM signal 
lives long enough to travel once through LISA's propa- 
gation arms places a constraint on the mass range of the 
SMBBH candidates: a few lO 5 A/ to (10 8 - 1O 9 )M . 

A step beyond detection analysis is that one of pa- 
rameter estimation. In Ref. [22|, [24[ it was found that 
through a single QNM detection LISA would be able to 
accurately infer the mass and spin of supermassive black 
holes: for black holes with mass M > 10 5 Mq the errors in 
mass and spin would be smaller than one part in 10 2 , and 
smaller than one part in 10 5 for the more optimistic case 
M > 5 x 10 5 A/ Q 2 . These predicted accuracies depend 
on the ringdown efficiency e r d, defined as the fraction 
of mass radiated in ringdown waves. In these references 
very conservative values were used: e r j ~ 0.1% — 3%. For 
example, it has been found in numerical simulations of 
two equal-mass, non-spinning black holes starting from 
quasi-circular motion that around e r( i ~ 2% — 3% of the 
total mass is radiated in the ringdown regime [25| . The 
inclusion of different masses and/or spin increases this 
value (see, for example, [26^. 

In [22|, [24[ it was also found that at least a second de- 
tection of either mass or spin should be possible for LISA. 
Resolving both spin and angular momentum (or, equiva- 
lcntly, both frequency and damping times associated with 
the QNM oscillation of this second mode) might require 
a very large critical signal to noise ratio, which might in 
turn need the second mode to radiate a significant por- 
tion of the emitted gravitational wave when compared to 
the first one. Whether this is feasible or not can only 
be established by giving precise predictions of the ampli- 
tudes for secondary candidates. 

Underlying in all these analyses is the implicit assump- 
tion that quasinormal modes and their spectrum of asso- 
ciated frequencies accurately describe the (intermediate) 
stage of the ringdown to a final Kerr stationary state. 
Which is certainly the case in linearized theory, but at the 
same time there is evidence that effects of self-interaction 
in gravitational waves might be observable with the ex- 
pected sensitivity of LISA [27j ■ 

Similarly, quasinormal modes also play an important 
role in the semi-analytical modeling of intermediate mass 
black holes (IMBH), such as in the Effective One Body 
approach (28l . |29| , where the gravitational wave as mod- 
eled within this formalism is stitched to a ringdown one 
consisting of QNMs by enforcing continuity of the wave 
and its derivatives and using the values of quasinormal 
frequencies as expected from linearized theory (30l - [32l | . In 
this context, it is worth recalling that in close-limit stud- 



Only cases with M > 1O 5 M are considered in [22I , I22I , |24H be- 
cause otherwise the QNM signal would be short lived enough 
that special detection techniques might be needed. 



ies of binary black holes it was found that corrections 
from second-order perturbations were in some cases sig- 
nificant [33|. For IMBHs, which could have total masses 
in the range of ~ lOOAf© — 10 Af , it is especially impor- 
tant to accurately model the merger and ringdown since 
they should fall in the frequency band of earth based 
gravitational wave detectors. Although the existence of 
IMBHs is still debatable, they could provide an interest- 
ing source for Advanced LIGO and VIRGO if they are 
present in dense globular clusters [Hj (see also ) . Re- 
cent observational evidence of an IMBH can be found in 

The previous discussions motivate us to study how 
nonlinear perturbations of black holes decay in time: do 
they do so just as in linearized theory or with a different 
spectrum of frequencies? We carry out our study through 
numerical simulations of first and second-order gauge- 
invariant gravitational perturbations of Schwarzschild 
black holes. We find that for all practical purposes, 
and to high numerical accuracy, the complex decay fre- 
quencies of second-order perturbations are the standard 
quasinormal ones from linearized theory -in contrast to 
previous predictions in the literature [40j-|42[~ and we ex- 
plain why this appears to be the case. Essentially, in all 
our simulations we find that mode-mode couplings ex- 
cite nonlinearities in the early stages of the perturbations 
before the quasinormal regime for the linearized pertur- 
bations kicks in. By the time the latter happens, those 
couplings have decreased to negligible values and the ex- 
cited nonlinearities essentially propagate as in linearized 
theory; and, in particular, they decay with the standard 
QNM frequencies. 

The structure of the paper is as follows. Section |TT] re- 
views the basics of Regge- Wheeler- Zerilli equations and 
quasinormal modes, and Sec. IHII the main features of the 
gauge- invariant approach here used for second-order per- 
turbations. Section ITVl describes our numerical approach 
and setup for solving the first and second-order Regge- 
Whcclcr and Zerilli equations, and Sec. [V] presents and 
discusses our main results. 



II. FIRST-ORDER PERTURBATIONS OF 
SCHWARZSCHILD AND QUASINORMAL 
MODES 

Metric gravitational perturbations can be expanded in 
tensor spherical harmonics. At the linear level modes 
with different angular structure decouple from each other 
if the background spacetime is spherically symmetric, as 
is the case for the Schwarzschild metric. In addition, 
for each multipole, perturbations of Schwarzschild can 
be further split into two sectors with different parities, 
which also decouple from each other in linearized theory. 

For each multipole, each of these sectors is completely 
described by a master function which depends on time 
and radius. The Regge- Wheeler function contains all 
the relevant information of the axial or odd-parity sec- 
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tor, whereas the Zerilli one encodes the polar or even 
degrees of freedom. These functions satisfy master evo- 
lution equations, 

□ ^>$7* -v RW il} $T = 0, (i) 
□ - y z {1 >*t* = o, (2) 

where {1 >v]/™ and {1} <I>™ denote, respectively, the first- 
order Zerilli and Regge- Wheeler functions for a given 
(I, m) mode. The potentials that appear in these equa- 
tions are given by 

_ 1(1+1) 6M r 2 A(A + 2) + 3A/(r - M) 
Z ~ ^ H (rA + Ml) 2 ' ^ ' 



and are, in particular, independent of the azimuthal 
index m. In these expressions A = \(l — + 2), 
M is the mass of the Schwarzschild black hole back- 
ground, r its areal radius, and □ is the two-dimensional 
D'Alambertian operator corresponding to the time and 
radial sector of the background. 

The complete spectrum of QNMs can be numerically 
obtained by analyzing Eqs. (|1I2|) in the frequency do- 
main. Computing in this way the amplitudes of QNMs 
for any given initial data is not so straightforward, 
though. Bearing in mind our motivation of studying the 
behavior of second-order perturbations, we instead solve 
Eqs. (|l|2j) in the time-domain. That is, we prescribe (a 
variety of) initial data, evolve them and analyze the solu- 
tions at different observer locations as functions of time. 

The early behavior of the solution depends on the type 
of initial data, followed by the QNM ringdown of the 
black hole. In Fig. [1] we show a typical solution of the 
Zerilli equation, for a fixed observer at r = 51.8M as a 
function of time. The initial data for this particular case 
consist of a Gaussian profile to the initial time deriva- 
tive of the Zerilli function, centered at r = 20M with a 
width er = AM, and the initial value of the Zerilli func- 
tion itself is set to zero. [This corresponds to what we call 
time-derivative initial data in the following sections, see 
Eqs. (|12I13|) ]. To measure the complex QNM frequency 
we perform a numerical fit to a function of the form 

f(t)=Ae at sm[b(t-t )} , (5) 

where the parameters to be fitted are A, a, b and to and 
the choice of the starting time for the QNM regime is 
chosen as that one which optimizes the fit, as introduced 
and explained in Ref. (43j . The quasinormal frequency 
is therefore given as w = a + hi. The expected value 
for the fundamental QNM for an I = 2 perturbation is 
0.37367 - 0.08896i [l5[, while our fit for this simulation 
yields 0.37077 - 0.08826i. 
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Figure 1: Solution to the first-order Zerilli equation at ob- 
server location r — 51. 8M. 



III. SECOND-ORDER GAUGE-INVARIANT 
PERTURBATIONS OF SCHWARZSCHILD 

We study second-order gravitational perturbations of 
Schwarzschild black holes using a gauge-invariant for- 
malism for arbitrary first and second-order perturbations 
[jjj . The key feature of this formalism is being able to 
consider perturbations with arbitrary angular multipolc 
structure, and has been possible mostly due to the devel- 
opment of a suitable theoretical framework [45|, H|| and 
to the advance of very efficient symbolic algebra tools for 
tensor- type calculations [13, EH • 

Next we very briefly summarize those results of the 
formalism presented in [44| which are relevant for the 
current work; see that reference for more details. 

Due to the intrinsic nonlincarities of General Rela- 
tivity, any non-trivial solutions of Eqs. (|1I2[) generate 
second-order contributions which are solutions of Zerilli 
and Regge- Wheeler-type equations with source terms, 

□ {2> $"' - ^rw l2> <f>? = {2} 5* , (6) 
□ f 2 %™ - v z {2} *r = {2} <S* . (7) 

The sources {2} <S* and {2} <S$ depend quadratically on 
the lower-order perturbations and their time and space 
derivatives from both first-order sectors. That is, in gen- 
eral the coupling of even (odd) parity modes generates 
odd (even) parity second-order modes; see Appendix lAl 
for a detailed description of the selection rules for this 
mode coupling. 

Second-order Regge- Wheeler- Zerilli (RWZ) type func- 
tions are not unique: one can add to them any quadratic 
combination of the first-order ones and they will still be 
gauge invariant and will still satisfy equations of the form 
(|6l7p with different sources. For this reason, asking what 
are the ringdown frequencies of second-order RWZ func- 
tions is not an unambiguous question; as opposed to, say, 
asking what are the ringdown frequencies of the emitted 
gravitational power. It turns out, however, that, as dis- 



4 



cussed in [44J and made explicit below, since the second- 
order quantities, like f2} *, {2} $, {2 >S* and {2} 5$, that 
we use correspond to the ones labeled as regularized in 
Ref. [HI, there is a simple relationship between emit- 
ted power and the RWZ functions. Under such choice of 
second-order RWZ functions it is therefore unambiguous 
and physically motivated to ask what are their ringdown 
frequencies. 

Reference [H| deals with the most general case for 
these sources and computation of the radiated energy. 
Here we quantitatively explore the ringdown frequencies 



for some particular cases. Namely, we study first-order 
{£ = \m\ = 2) even-parity and {£ = 2, m = 0) odd-parity 
perturbations and the modes that they generate through 
self-coupling. That is, for simplicity we do not consider 
the coupling between the mentioned even and odd-parity 
first-order modes. We do not explicitly list the sources 
{2} <S^ and {2} S<i, for the cases here considered because 
they are rather lengthy and complicated expressions, but 
they are available from the authors upon request. The 
generated modes and the radiated energy carried by them 
are described next and summarized in Table Q] 
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Table I: First-order modes considered and second-order ones generated by self- coupling. 



A. CASE A: even-parity I = \m\ = 2 perturbations 
and generated modes 

As discussed in Appendix [^1 the self-coupling between 
these modes generates second-order {£ = 4, m = ±4) 
even-parity (polar) ones, whereas the coupling between 
them (m = 2 with m = —2) gives rise to second-order 
(£ = 4,m = 0), (£ = 2,m = 0) and {I = 0,m = 0) 
even-parity (polar) modes, and {£ = 3, m = 0) and {£ = 
l,m = 0) odd-parity (axial) ones. Since this paper only 
deals with different radiative aspects of the system, we 
can ignore modes with £ < 2. 

Furthermore, we assume that the Zerilli functions 
{1} *± 2 that describe the first-order (I = 2,m = ±2) 
modes are real 



'"^el. (8) 

This means that both modes are described by the 



same function, since generically the relation ( {1} vl>™)* = 
(— l) m { 1 >\lrr' m holds, and, in essence, the system reduces 
to a problem of self-coupling. The second-order even- 
parity (polar) modes inherit this property, in such a way 
that we only need three functions to describe them: 

Due to the assumption ([5]), none of the second-order 
odd-parity (axial) modes are generated. This happens 
because the source for a m = mode must be real. 
Schematically, the generic term of the source for this ax- 
ial (I = 3,m = 0) mode can be written as {1} *^ 2 
and it is straightforward to see that its real part vanishes 
under the assumption (|8j). 

In this particular CASE A, the radiated power asso- 
ciated with the mentioned modes for a given observer 
located at r b s as a function of time is given by 



J 



Power(r 6 s ,t) 



dE 
~dt 
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127T 
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640 



|cV 2 >^| 2 



d t {2 m 2 } + ^\d t ^° 2 \ 2 + <D(e 5 
) yo7r 



(9) 



where e is the pcrturbativc parameter and all the expres- 
sions on the right-hand side are evaluated, naturally, at 
( r obs,t). In principle this equation is valid only at null 
infinity but, as it is usually the case in computations, we 
evaluate it at a finite radius. 



B. CASE B: odd-parity (£ = 2, m = 0) perturbations 
and generated modes 

In principle, following the selection rules, the odd par- 
ity (£ = 2, 77i = 0) mode would generate by self-coupling 
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the second-order {£ = 0, m — 0), [l = 2, m = 0) 
and (£ = 4, m = 0) even-parity modes, as well as the 
(I = 1, m = 0) and (i = 3, m = 0) odd-parity ones. How- 
ever, in this axisymmetric case (m = for all modes), 
the Clebsch-Gordan-like coefficients that appear in the 
sources {21 5$ and {2> <S* vanish for those second-order 
modes with an odd harmonic label I (see Appendix fA| . 



IV. NUMERICAL APPROACH FOR SOLVING 
THE MASTER EQUATIONS 

We now describe in some detail our numerical approach 
for evolving the first and second-order RWZ equations, 
since in the past difficulties have been reported with the 
high-order derivatives in the sources of these equations. 

We numerically solve the first and second-order equa- 
tions using a pseudo-spectral collocation method. The 
spatial derivatives are computed using Chebyshev poly- 
nomials and Gauss-Lobatto (GL) collocation points, and 
the system is evolved in time using a standard fourth- 
order Runge-Kutta scheme. We use a small enough 
timestep for the time integration so that the solution 
converges exponentially with the number of collocation 
points (see below). The accuracy of all the simulations 
presented in this paper are at the level of double precision 
round-off. 

GL collocation points are not equally spaced; instead, 
they cluster near the edges of the computational domain 
(equally spaced points would not give exponential con- 
vergence). For that reason it is standard to use a multi- 
domain approach. Here we subdivide our radial domain 
in (non-overlapping) blocks of length 10M each, commu- 
nicated through a penalty technique. At the interface 
each incoming characteristic mode u + is penalized ac- 
cording to (see [49| and references therein) 

• aNH + + 
u + = (...) [u — v ) 

^block 

where v + is the value of the same mode at the interface 
point using the neighboring block, r^i oc ^ is the size of 
the corresponding block (10M in these simulations), a 
is the associated characteristic speed, TV the number of 
collocation points on that block and 8 a penalty param- 
eter chosen here to be 5 = 0.6. At the outer boundary 
each characteristic incoming mode is similarly penalized 
to zero; though this is done simply to achieve stability, in 
our simulations the domain is large enough that our re- 
sults are causally disconnected from the outer boundary. 
The singularity of the black hole is dealt with through 
excision. That is, by using Kerr-Schild coordinates for 



Hence, none of the second-order odd-parity modes are ex- 
cited and, up to this order, we are left with two radiative 
second-order modes, 

{{ 2 >$°, f 2 >*°}. 
The radiated power is then given by 



(10) 

I 

the background spacetime and placing an inner bound- 
ary inside the event horizon. 

As an example, Fig.[2]shows a self-convergence test for 
the first-order {1 = 2 = in) Zcrilli function, extracted at 
r = 5 1.8 M 3 , both changing the number of collocation 
points as well as the timestep. The initial data used for 
such test correspond to the same one used in Fig. [TJ 

{1)^2 = ; = ,r) = e~ {r - ro)2/a2 , (11) 

with a = AM, rg = 20M, and the spatial domain r G 
[1.8M,301.8M]. In Fig. [3] we also show the result of a 
convergence test for the generated second order (I = 4 = 
m) Zerilli mode. From these figures we see that using 
30 collocation points per domain and a timestep At = 
0.001M gives a numerical error at the level of double 
precision round-off; from hereon we use such resolutions 
for all our simulations. 

In order to compare the magnitude of the errors with 
the solutions themselves, in Fig. 2] we show the ab- 
solute values of the first-order {1> \1/| and second-order 
|{2}^ ( ) {2K[/4j. Zerilli solutions from the previous 

plots at their highest resolutions, all extracted at the 
same observer location. 

We note in passing that for most of the ringdown the 
order of magnitude of the second-order Zerilli functions 
appears to be comparable to (and in one case even larger 
than) the first-order one. There is no contradiction in 
this, since their contribution to the radiated energy is 
scaled by e 4 , while the contribution of the first order Zer- 
illi function is scaled by e 2 , sec Eq. ©. 

A. Setup of numerical simulations 

We could introduce non- vanishing second-order modes 
via initial second-order perturbations. However, we are 
interested in mode-mode coupling. Put differently, we 



3 We place observers at the beginning/end of each domain: 
1.8M, 11.8M, 21.8M, etc. 
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Figure 2: Numerical errors for different spatial resolutions 
using a fixed timestep At = O.OfM (top), and for different 
timesteps using a fixed spatial resolution of TV = 60 points per 
domain (bottom). Both figures show the differences between 
several resolutions and the most accurate one, which is TV — 
60 for the top panel and At4 = 0.0005M for the bottom one. 
In both cases the observer is located at r = 51.8M. We 
see exponential convergence and errors in the order of double 
precision round-off. 



Figure 3: Numerical errors in the second-order Zerilli function 
{2} ^| for different spatial resolutions and a fixed timestep 
At — 0.001M. The errors are to be interpreted as in the 
previous figures. For TV — 30 (which are the ones used in the 
rest of this paper) and higher, they are at the level of double 
precision round-off. 
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Figure 4: Absolute value of first and second-order Zerilli func- 
tions extracted at r = 51.8M. 



are interested in the particular solution of Eqs. (|7I6[) (that 
is, the one with vanishing initial data), since the homo- 
geneous one will be exactly the same as at first order. 
Therefore, in this paper wc always impose vanishing ini- 
tial data for all the second-order modes and concentrate 
on those modes generated by first-order mode coupling. 

In the following section we solve the first-order RWZ 
equations with four different types of initial data, varying 
both the location ro and width a of the initial data, 

1. Time Derivative (TD) 

{1 >*2(t = 0,r) = 0, (12) 
{1} *2( t = ,r) = e-t—o) 2 /- 2 . ( 13 ) 



2. Time Symmetric (TS) 

fl, *2(t = ,r) = Me- {r - ro ^ /<T \ 
^l(t = 0,r) = 0. 

3. Approximately Outgoing (OUT) 

{11 *2(t = ,r) = Me-( r - r »» 2 / ff2 , 

^H>l(t = 0,r) = -(l-2M/r)d r {1} ^>l(t = 0,r). 

4. Approximately Ingoing (IN) 

fl} ^(i = 0,r) = Me- {r ~ ro)2/a2 , 

{1 %2( t = ,r) = (l-2M/r)a r 11| ^(t = 0,r). 
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V. DECAY FREQUENCIES OF 
SECOND-ORDER PERTURBATIONS 

Ioka and Nakano have put forward the suggestion that 
at second order new frequencies should appear in the 
ringdown spectrum, which would be given by the sum of 
different pairs of standard QNM frequencies. In particu- 
lar, according to this prediction, the dominant frequency 
would correspond to the double of the standard funda- 
mental one from linearized theory [|o|, |4l| . This seems 
reasonable, since the sources for the second-order mas- 
ter equations are quadratic in the first-order modes, so 
one might well expect that frequencies get summed-up 
in Fourier space. The physical picture, however, appears 
to be at the same time more subtle and simpler: our nu- 
merical simulations indicate that in practice second-order 
perturbations decay with the standard QNM frequencies 
from linearized theory. 

Recall that the physical process here studied is the cou- 
pling of linear modes. That is, we initialize all second- 
order perturbations to zero. The second-order mas- 
ter equations have sources which are, indeed, quadratic 
in the first-order modes. What we observe in all our 
simulations, though, is that those sources quickly ex- 
cite the second-order perturbations and afterwards decay 
very fast in time. As a consequence, once the second- 
order modes have been excited and reached the regime 
in which they oscillate with a constant complex fre- 
quency (i.e. what in linearized theory corresponds to the 
QNM regime), they essentially propagate with a vanish- 
ing source. In other words, as first-order perturbations 
do. And, in particular, they do oscillate and decay with 
the same, standard, QNM frequencies from linearized 
perturbation theory. 

We show this behavior in some detail for the four ini- 
tial data types (Sec. IIV A[) of CASE A perturbations 
(Sec. lIII A|) . with r = 20M and a = AM. Figure [5] shows 
the first-order Zerilli function and the (£ = 2,m = 0) 
second-order one for different initial data profiles. In all 
cases the source decays much faster than the second-order 
solution itself and therefore its role in determining the be- 
havior of the latter by the time it enters the QNM regime 
is negligible. From the same figure one can notice that 
the type of initial perturbation does determine the time 
by which the second-order Zerilli function enters the tail 
regime; but this is not surprising, since it already hap- 
pens for the first-order one. 

In order to gain further insight into these observations 
we display the dynamics of first and second-order Zerilli 
functions and the source term {2} Sy of Eq. (J7J, now for 
the [l = A, m = 0) second-order mode. Fig. [6] shows 
these three quantities as functions of radius at different 
times. The source term is dominant only during the first 
~ 20M, later decaying at a fast rate to several orders of 
magnitude below the second-order Zerilli function. 

Table |TT] shows the fitted QNM frequencies from our 
numerical data, for the different initial-data profiles of 
CASE A perturbations, with r = 20M and a = AM. 
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Figure 5: CASE A simulations: first-order (£ = 2, m = 2) 
Zerilli function and second-order (I = 2, m = 0) one, along 
with the source term for the second-order master equation 
[Eq. l|7]l]. as functions of time for different initial data profiles. 
The source plays a role only at very early times in exciting 
the second-order modes, afterwards being much smaller than 
the first and second- order solutions. 
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Figure 6: CASE A simulations: snapshots of the first-order 
Zerilli function and second order (£ = 4, m = 0) one, along 
with the source term, for ingoing initial data. The generic 
behavior of the source is to rapidly decrease several orders of 
magnitude below the solutions themselves. [Note: In the first 
snapshot the second-order Zerilli function vanishes because it 
corresponds to t = and our setup of the physical problem.] 



They agree quite well with those predicted by first-order 
theory for each of those modes. 

As described in Sec. IIVI our numerical simulations are 
of very high accuracy (both at first and second-order) : all 
the errors are at the level of double precision round-off, 
~ Kr 14 -1CT 12 (see Figs.03]). Therefore, we do not con- 
sider lack of resolution as a possible reason for not finding 
traces of the predicted new second-order QNM frequcn- 



lst order 2nd order 



ID £ = 2, m = 2 err % 


£ = 2, m = err % 


TD 0.37077 - 0.08826; 0.8 
TS 0.37353 - 0.08837; 0.2 
IN 0.37061 - 0.08887; 0.8 
OUT 0.37107 - 0.08624z 1.0 


0.37334 - 0.08883i 0.1 
0.37335 - 0.08766; 0.3 
0.37373 - 0.08945; 0.1 
0.37074 - 0.08902; 0.8 




2nd 

ID £ = 4, m = err % 


Drder 

£ = 4, m = 4 err % 


TD 0.80916 - 0.09418i 0.003 
TS 0.80920 - 0.09420; 0.005 
IN 0.80918 - 0.09416; 
OUT 0.80931 - 0.09425i 0.019 


0.80916 - 0.09418; 0.003 
0.80920 - 0.09420; 0.005 
0.80918 - 0.09416; 
0.80931 - 0.09425; 0.019 



Table II: Measured quasinormal frequencies from our numeri- 
cal simulations (CASE A). They agree with those predicted by 
linearized theory, even for the second-order modes generated 
due to mode-mode coupling. The predicted QNM frequencies 
from standard linearized perturbation theory, as quoted in 
QJ], are 0.37367 - 0.08896; for £ = 2 and 0.80918 - 0.09416; 
for £ = 4 (in linearized theory there is degeneracy with re- 
spect to the azimuthal index m. The relative errors for real 
and imaginary parts of the measured frequencies u are com- 
puted as |w-w e xactl/l w exactl- 



cies in our simulations. Similarly, one might think that 
those predicted frequencies are in fact present, but with a 
very small amplitude. Figure [7] indicates that in practice 
that does not seem to be the case: the residual of the fit 
for the second-order Zerilli function [in the case shown it 
is the [I = 2, m = 0) one, for TD initial data, first-order 
perturbations] — defined as the function minus its fit — 
does not appear at all to correspond to an oscillation 
and decay with twice the standard complex fundamental 
quasinormal frequency for that mode (or to an overtone). 
Instead, it appears to be the residual associated with the 
fact that QNMs are not complete. 

For completeness, in Appendix [B] we provide results 
of the fitted frequencies for the four initial data types of 
CASE A perturbations, now varying both the location 
and width of the initial data; all of them support the 
same conclusion. 

Finally, we briefly discuss the results of some CASE 
B [odd-parity (I = 2, m = 0) first-order mode] pertur- 
bations, since the conclusions are identical. As discussed 
in Sec. IIII Bl and summarized in Table HI they generate 
both {I = 4, m = 0) and (£ = 2, m = 0) even-parity 
second-order modes. The fitted frequencies from the nu- 
merical solutions (for a simulation of TD linear initial 
data with r$ = 20M and a = 4M) for these second- 
order modes yield, respectively, 0.37441 — 0.08921i and 
0.80932 — 0.09419i, to be compared with the expected 
values from perturbation theory: 0.37367 — 0.08896; for 
£ = 2 and 0.80918 - 0.09416; for I = 4. 
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Figure 7: Residual from fitting a second-order Zerilli function 
to a complex frequency mode. The fitted frequency corre- 
sponds to the standard fundamental quasinormal one for that 
multipole index, and the residual does not appear to contain 
traces of twice that frequency (or any other). 

VI. FINAL REMARKS 

In this paper we have numerically evolved first and 
second-order self-generated gauge-invariant gravitational 
perturbations of Schwarzschild black holes with a variety 
of initial data sets, studying the oscillation and decay be- 
havior of nonlinear modes and, more specifically, whether 
they correspond to the standard QNM frequencies or to 
a different spectrum. We have found, in all cases, that 
second-order modes decay through the standard QNM 
frequencies, and that the picture behind this is remark- 
ably simple: first-order perturbations trigger high-order 
ones through source terms which afterwards rapidly de- 
cay in time. Besides, by the time the solutions reach the 
regime in which they oscillate and decay at a constant 
rate (the QNM regime in the case of linearized perturba- 
tions), the second-order modes for all practical purposes 
propagate as in linearized theory. 

Mode-mode coupling in the ringdown of black holes has 
been previously studied through numerical simulations of 
full Einstein equations [HO, HH; however, no conclusions 
seem to have reached or sought for in terms of the devi- 
ations in the ringdown spectrum from the linearized one 
(presumably due to lack of resolution). 

The fact that nonlinear aspects of Einstein's equations 
in the ringdown of black holes appear to already be cap- 
tured — at least in what oscillation and decay frequencies 
concerns — by linearized perturbation theory is somehow 
remarkable and should be of use both for modeling black 
holes in the ringdown regime as well as in data analysis 
searches of gravitational waves. 
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Appendix A: Mode coupling: selection rules 

In this appendix we summarize, for completeness, the 
selection rules for mode coupling (4f||46|. A second-order 
(I, m)-mode gets a contribution from a pair of first-order 
modes (l,rh) and (I, in) if three conditions are obeyed. 
First, the harmonic labels must be related by the usual 
composition formulas 

\l - l\ < l < l + l, and m + m = m. (Al) 



Second, mode coupling must conserve parity. To any har- 
monic coefficient with label I, we associate a polarity sign 
a such that, under parity, the harmonic changes by a sign 
er(— 1) . Polar/even parity (axial/odd parity) harmonics 
have a = +1 (a = —1). Then, parity conservation im- 
plies the third condition: 



(-l) T+l ~ l = aaa, (A2) 



where a and a are the polarity signs corresponding to 
the modes (l,m) and (Z,m) respectively. In particular, 
there is a special case in which the coupling of two modes 
satisfying Eqs. (|A1[) and (|A2I) does not contribute to a 
second-order mode, and the reason comes from the prop- 
erties of the Clebsch-Gordan-like coefficients that ap pea r 
in the product formula for the tensor harmonics [45j. 
In axisymmetry (m = m = 0) the mentioned Clebsch- 
Gordan-likc coefficients vanish if I + I + I is odd. 

This analysis can be extended to higher orders. In 
particular, the parity condition implies that a collection 
of k modes with harmonic labels {l±, and polar- 
ities {o"i, (Tfc} contributes to the mode (l,u) only if 

(-i)V = nt 1 (-i)' 4 ^. 
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Appendix B: Numerical decay frequencies of second-order perturbations from time domain simulations 







ro £ = 2, m 

20 0.37373 

40 0.37147 

60 0.37141 

80 0.37102 

100 0.37101 - 0.08840; 



Ingoing 
I = 4, m = 
0.08949i 0.80918 - 0.09416; 
0.08833i 0.80916 - 0.09417; 

0.09413; 



0.08834; 0.80915 
0.08839; 0.80917 



£ = 4, m = 4 
0.80918 - 0.09416; 
0.80916 - 0.09417; 
0.80915 - 0.09413; 



0.09415; 0.80917 - 0.09415; 
0.80888 - 0.09414; 0.80872 - 0.09417i 



£ = 2, m = 
0.37074 - 0.08902; 
0.37079 - 0.08891; 
0.37093 - 0.08844; 



Outgoing 
£ = 4, m = 
0.80931 - 0.09425; 
0.80911 - 0.09409; 
0.80911 - 0.09417; 



0.37140 - 0.08816; 0.80913 - 0.09413; 
0.37148 - 0.08782; 0.80914 - 0.09426; 



£ = 4, m = 4 
0.80931 - 0.09425; 
0.80911 - 0.09409; 
0.80911 - 0.09417; 
0.80913 - 0.09413; 
0.80914 - 0.09426; 



a 

2 0.37708 - 0.09116; 

4 0.37373 - 0.08949; 

8 0.37286 - 0.08893; 



0.80918 - 0.09418; 
0.80919 - 0.09417; 
0.80875 - 0.09407; 



0.80918 - 0.09418; 
0.80918 - 0.09417; 
0.80875 - 0.09407; 



0.37365 - 0.08891; 
0.37079 - 0.08902; 
0.36887 - 0.08284; 



0.80918 - 0.09416; 
0.80931 - 0.09425; 
0.80923 - 0.09599; 



0.80918 - 0.09416; 
0.80911 - 0.09409; 
0.80921 - 0.09599; 



r t = 2, m = 

20 0.37334 - 0.08883; 

40 0.37277 - 0.08893; 

60 0.37378 - 0.08921; 

80 0.37383 - 0.08995; 

100 0.37387 - 0.08999; 



Time Derivative 
£ = 4, m = 
0.80916 - 0.09418; 
0.80919 - 0.09417; 
0.80920 - 0.09418; 
0.80915 - 0.09411; 
0.80914 - 0.09409; 



£ = 4, m = 4 
0.80916 - 0.09418; 
0.80919 - 0.09417; 
0.80920 - 0.09418; 
0.80915 - 0.09411; 
0.80914 - 0.09409; 



£ = 2, m = 
0.37335 - 0.08766; 
0.37376 - 0.08984; 
0.37267 - 0.08948; 
0.37344 - 0.08852; 
0.37321 - 0.08823; 



Time symmetric 
£ = 4, m = 
0.80920 - 0.09420; 
0.80915 - 0.09413; 
0.80917 - 0.09416; 
0.80917 - 0.09413; 
0.80923 - 0.09422; 



£ = 4, m = 4 
0.80920 - 0.09420; 
0.80915 - 0.09413; 
0.80917 - 0.09416; 
0.80917 - 0.09413; 
0.80923 - 0.09422; 



2 0.37360 - 0.08906; 0.80917 - 0.09417; 0.80917 - 0.09417; 
4 0.37334 - 0.08893; 0.80916 - 0.09417; 0.80916 - 0.09418; 
8 0.37315 - 0.08935; 0.80924 - 0.09342; 0.80924 - 0.09334; 



0.37404 - 0.08972; 0.80917 - 0.09416; 0.80916 - 0.09416; 
0.37335 - 0.08766; 0.80920 - 0.09420; 0.80920 - 0.09420; 
0.37943 - 0.08760; 0.80867 - 0.09424; 0.80867 - 0.09423; 



Table III: For £ = 2 the expected fundamental frequency from linearized theory is 0.37367 — 0.08896;, while for £ = 4 frequency 
it is 0.80918 - 0.09416; hjl. Shown are the fitted values from our numerical solutions to the second order Zerilli equations, 
for four different sets of initial data. For each initial data configuration we vary the location ro and width a of the linearized 
perturbations, keeping one fixed (ro = 20M and a = 4M, respectively) while varying the other. The fitted frequencies appear 
to be insensitive to the choice of initial data and agree with the QNM frequencies from linearized perturbation theory. 
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